Nonlinear intraband tunneling of BEC in a cubic three-dimensional lattice 



V. S. Shchesnovich^ and V. V. Konotop^ 
^ Instituto de Fisica - Universidade Federal de Alagoas, Maceio AL 57072-970, Brazil 
^ Centra de Fisica Teorica e Computacional, Universidade de Lisboa, 
Complexo Interdisciplinar, Avenida Professor Gama Pinto 2, Lisboa 1649-003, 
Portugal; Departamento de Fisica, Faculdade de Ciencias, Universidade de Lisboa, 
Campo Grande, Ed. C8, Piso 6, Lisboa 1749-016, Portugal 

The intra-band tunneling of a Bose-Einstein condensate between three degenerate high-symmetry 
X-points of the Brillouin zone of a cubic optical lattice is studied in the quantum regime by reduction 
to a three-mode model. The mean-field approximation of the deduced model is described. Compared 
to the previously reported two-dimensional (2D) case [Phys. Rev. A 75, 063628 (2007)], which is 
reducible to the two-mode model, in the case under consideration there exist a number of new stable 
stationary atomic distributions between the X-points and a new critical lattice parameter. The 
quantum collapses and revivals of the atomic population dynamics are absent for the experimentally 
realizable time span. The 2D stationary configurations, embedded into the 3D lattice, turn out to 
be always unstable, while existence of a stable ID distribution, where all atoms populate only one 
X-state, may serve as a starting point in the experimental study of the nonlinear tunneling in the 
3D lattice. 

PACS numbers: 03.75.Lm 



I. INTRODUCTION 

In recent paper |T] it has been shown that exploring 
the nonlinear tunneling of a Bose-Einstein condensate 
(BEC) loaded in a two-dimensional (2D) optical lattice 
allows for theoretical and experimental study of diversity 
of fundamental issues of the nonlinear physics. Among 
them we mention the validity of the mean-field approx- 
imation (considered previously in Ref. [2]), which was 
used in the study of the nonlinear tunneling [3J H], the 
accuracy of the semi-classical (i.e. WKB) approxima- 
tion, where the inverse number of atoms plays the 
role of the effective Planck constant (similar to the ideas 
reported in Ref. [5j for a coupled two-mode model and 
in Ref. [B] for the Bohr-Sommerfeld quantization rule), 
and the macroscopic manifestation of the quantum col- 
lapses and revivals - a pure quantum effect related to the 
discrete nature of quantum spectra [7] . The approach de- 
veloped in Ref. [1^ was based on the energy degeneracy 
due to the rotational, more specifically C4, symmetry of 
the lattice, where the modulationally stable Bloch states 
at the high-symmetry A"-points were used (for either at- 
tractive BEC loaded in the first lowest band or repul- 
sive BEC loaded in the second band). The modulational 
stability allows one to reduce the mean-field description 
of spatially inhomogeneous matter waves to the effective 
two-mode model describing the populations of the reso- 
nant states (see also [8 and the references therein). It 
was shown in Ref. HJ by direct numerical simulations that 
the two-mode model gives a remarkably well description 
of the dynamics for a relatively long interval of time. 

Since optical lattices are routinely available in all di- 
mensions (see, for instance, the recent review [9J), an 
intriguing possibility is to explore the dynamics similar 
that reported in Ref. [1] but in the 3D case. Consider- 
ing a cubic lattice one expects the dynamics to be much 



richer as compared to the 2D case, because the distinct 
X-points of the same band are now three-fold degener- 
ate. Concentrating on the modulationally stable case, 
the only situation considered in the present paper, one 
would expect that the stable distributions are very dif- 
ferent from those in the 2D configuration, moreover, the 
latter (embedded in the 3D lattice) is found to be un- 
stable. The most significant feature of the problem at 
hand is that it reduces to an effective three-mode model, 
whose dynamics is described by a Hamiltonian with two 
degrees of freedom, due to the constraint imposed by con- 
servation of the number of atoms. Taking into account 
that such dynamics in the vicinity of a stable point is, 
generally speaking, characterized by two frequencies, as 
well as the facts that in most of the spectrum the energy 
level spacing in the quantum model scales at least as iV~^ 
(except, for example, the local bound states close to the 
semiclassical stationary points) and that the number of 
atoms used in BEC experiments is large, one can expect 
that the phenomenon of quantum collapses and revivals 
in the respective quantum system is significantly affected 
(and even suppressed completely). At the same time new 
features can be expected due to the fact that the classical 
motion now can be either regular (in the vicinity of the 
stationary points) or chaotic, what will naturally affect 
the underlying quantum evolution. 

Study of the phenomena mentioned above with spe- 
cial attention payed to the correspondence between quan- 
tum and semi-classical dynamics (like in the 2D case, the 
semi-classical dynamics will be obtained in the mean-field 
approach) constitutes the main goal of the present paper. 

More specifically, we start by introducing in Sec. |TT] 
the quantum model and discussing its validity, physical 
parameters and the time span achievable in possible ex- 
periments. The mean-field limit is derived in Sec. |III| by 
making use of the WKB approximation for the quantum 
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model rewritten as a Schrodinger equation for an effec- 
tive quantum particle. We study the stationary points of 
the mean-field dynamics and their stability by locally di- 
agonalizing the Hamiltonian. In Sec. |IV| we compare the 
numerical simulations of the quantum model to those of 
the mean-field approach. Finally, in Sec.|V]we summarize 
the results. 



II. THE REDUCED QUANTUM MODEL 

A. The reduced Hamiltonian 

Let us start with the Hamiltonian of a BEC in an op- 
tical lattice 

J dW(x) + F(x)^ i/^(x) 

V 

+ ^ J dW(x)^^(x)^(x)V^(x), (1) 
V 

where x e M'^, V{x) is the cubic optical lattice potential, 
V = Mvq is the total volume of the lattice consisting of 
M cells each one of the volume vq, g is the interaction 
coefficient, and and tpi^) ^^'^ the creation and an- 

nihilation field operators. Introducing the Bloch waves 
<y5nk(x) through the standard eigenvalue problem 



2m 

where n is a number of the zone and k is the wave- vector 
in the first Brillouin zone (BZ), we expand 



(2) 



By analogy with the 2D case [T] , after the symmetry of 
the lattice is fixed the nonlinear tunneling phenomenon 
in 3D does not depend on the particular shape of the 
potential, as for instance on its separability: the only 
lattice parameter A entering the model, see Eq. (12) be- 



low, characterizes the effective lattice depth rather than 
its geometric properties. Therefore we concentrate on the 
simplest case of a separable cubic lattice: 



V — sEr[cos{x/(l) + cos(?//d) -t- cos(z/d)], 



(5) 



where d is the lattice period and the lattice depth s 
is measured in the units of the recoil 
h'^TT'^/{2m(f). The respective BZ is given by [-5, f ] x 

We consider the three distinct X-points, Xi = 
(5,0,0), X2 = (0,f,0) and X3 = (0,0, f), degenerate 
in the Bloch energy, which pertain to the modulationally 
stable Bloch band, and assume that only these points 
are significantly populated by the BEC atoms at t = 0. 
Then, repeating the arguments of the 2D case [T] it can 
be shown that for small g (see also below) the energy and 
quasi-momentum conservation laws allow one to discard 
the transitions, due to the scattering of BEC atoms, to 
all other points except the transitions between the three 
degenerate X-points. We arrive at the three-mode ap- 
proximation: 



V'(x) = V7l(x)6i + (^2(x)fe2 + </?3(x)&3, 



(6) 



where the Bloch function (/5j(x) and the operator bj cor- 
respond to the Xj-point (here we use the simplified-index 
notations for the populated states). Substitution of this 
expression into Eq. ([T]) results in the approximate Hamil- 
tonian 



Xiblb]b,b, 



where the creation and annihilation operators satisfy 
the usual commutation relations [fe„k, ^n'k'] = and 
[^Tik: &|i/k'] — ^rm"^kk'- The cxpausion ^ allows one to 
rewrite the Hamiltonian ([ij in the form 



n,k 

E nin2ri3n4 r hi^ h h 



ki,...,k4 



(3) 



where Q is an arbitrary vector of the reciprocal lattice 
and 

= f / d^X<^k^<^k^^„3k3^„,k. (4) 



(hereafter the asterisk stands for the complex conjuga- 
tion). 



x2\4Y.¥kb,bk+Y.(b]nb,r\ (?) 
j<k J 

where Ex is the Bloch energy at the X-point and the 
only nozero coefficients (El are given by 



Xi 
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jd^^WlW X2^\j d^^Wl?W2?- (8) 



The Hamiltonian ^ commutes with the total number of 
atoms N in the X-points 



(9) 



which reflects the approximation made. 

The details of the derivation can be found in Ref. [1]. 
Let us, however, present an alternative way to arrive at 
the Hamiltonian ([?]) . First we note that if initially only 
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the X-points of the same Bloch band are populated, the 
rates of the quantum transitions to other points (within 
the same band and to other Bloch bands) due to the s- 
wave atomic scattering (i.e. the nonlinearity of BEC), 
treated as a perturbation, are defined, according to the 
Fermi golden rule, by the energy conservation between 
the initial and final states within the linear model and 
are proportional to the density of the states at the par- 
ticular point to which the transition takes place. In our 
particular case, an additional rule on the transitions, sat- 
isfied by the nonlinearity of BEC, is that the sum of the 
quasi-momenta of the two atoms before the scattering 
and after that can differ only by some reciprocal lat- 
tice vector (we consider the Bloch waves as the basis). 
Recalling that the density of states is large only on the 
boundary of the Brillouin zone (it diverges in the limit 
of an infinite lattice) and only the X-points of the same 
Bloch band on the boundary have the same energy and 
satisfy the quasi-momentum conservation (in the above 
sense), all transitions except those to the X-points of 
the same band can be neglected. The same conclusion 
is also derived within the mean-field approach. Indeed, 
the resonant four-wave processes are determined by the 
phase-matching conditions (equivalent to the energy and 
the quasi-momentum conservation) and the population 
of the respective points, thus unpopulated points (or en- 
tire bands) which do not satisfy the matching conditions 
do not make any contribution to the resonant processes. 

In this set of arguments, the transitions to other (res- 
onant) points of the same Bloch band are neglected due 
to either the negligibly low density of states compared 
to that at the boundary points or due to the quasi- 
momentum non-conservation. On the other hand, the 
transitions to the boundary points of other Bloch bands 
are neglected due to the energy conservation, i.e. un- 
der the condition that the nonlinearity of BEC is much 
smaller than the band gap at the boundary of the Bril- 
louin zone. In the shallow lattice limit s <C 1 (A close 
to 1), for instance, the condition of a relatively large gap 
requires that Xi^ ^ 9^1^ ^ Er- 

On the other hand, we use the expansion over the 
Bloch-wave basis [see Eq. ([2|], what means that a few- 
mode approximation implies concentration of particles in 
the respective states. This imposes a constraint on the 
potential depth. Indeed, two body interactions result in 
nonzero spectral width of a Bloch state, which in our case 
can be estimated as ujnl — In order to be able to 

neglect the effect of the spectral width on the dynamics 
(what is done in the present paper) one has to require 
it to be much less than the spectral width of the lowest 
band (the most narrow band in the general situation). 
The latter is typically of order of the recoil energy Er- 
Thus we require watl Er/h (notice that that this is 
the limit in some sense opposite to the standard condi- 
tions of applicability of the Bose-Hubbard model, where 
due to relatively large amplitude of the periodic potential 
and, consequently, strong on-site localization of the wave 
functions, the expansion over Wannier functions is more 



appropriate). In order to get an idea about the physical 
range of the parameters, let us estimate the frequency 
ujnl- The coefficient Xi — 2 I '^^'^Ifl^ hence we 

can estimate 



7 



Snhas N _ 1 as N 
m V 7 m V ' 
3.78 X lO^^{J-sy\ 



(10) 



For instance, for ^^Rb we have m — 1.44 x and 
fls = b.lnm thus for a lattice with M ~ 20'^ cells with the 
lattice constant d — Ijim (wq = 10~^^cm'^) and N — 1000 
we get ojnl ~ Vl.bHz. On the other hand assuming the 
potential depth to be equal to the recoil energy we obtain 
the width of the first lowest band to be 0.59 • E^/h « 
590 Hz, and thus fkoNL/Er « 0.021. By reducing the 
potential depth this relation can be further improved. 

In this context it is also relevant to mention, that 
the dimensional time r is measured in the units I/lonl, 
which for the above example is approximately 0.08 s. 
Taking into account that the characteristic lifetime of 
a condensate can reach 10s, we conclude that for typical 
experiments the characteristic time is about 100 dimen- 
sionless units. This time can be significantly enlarged 
(i.e. making observation of the reported effect much eas- 
ier) by using lighter, say lithium, atoms and/or with a 
larger s-wave scattering length, achievable by Feshbach 
resonance. 



B. The dynamical equations 

It is convenient to use the dimensionless time r = 
uJNLt, subtract from the Hamiltonian ([t]) the constant 
term Hq — {Ex — Xi)^ [here we use Eq. ^] and nor- 
malize the result by dividing by N'^ , which allows for the 
transition to the mean- field limit — s- 00, since the re- 
sulting Hamiltonian is written in the population densities 
fij/N. This results in the Schrodinger equation 



TV 

with the Hamiltonian 



(11) 



H 



f4E".^+A5:n,n, + ^E(^t)vA(12) 



where A — XilX\- Equation (12 1 may be interpreted as 
a Schrodinger equation for a single quantum particle (see 
also Eq. ([^ below). 

Denoting by kj the number of atoms populating 



the Xj-point (such that fci 



N) one can 



expand the wave function l^*) over the Fock basis 
|fci,/c2,iV-fci-fc2) = |fci)|fc2)|A^-A:i -fe): 

E E C'fc„fe,(i)|fci,A:2,iV-fci-fc2), (13) 

fcl=0 fc2=0 
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where the expansion coefficients obey the normalization Now Eq. (12) can be cast in the form: 
condition 



N N-ki 
ki=0 k2=0 



(14) 



N dr 



where 



dki-l,k2 + lCki-2,k2+2 + dki + l,k2-lCki+2,k2-2) i (15) 



ak,M - 1 + 2(2A - l)A^-2 [(^^ ^ „ _ _^ ^ 

&fei,fe2 = {fci(l + ki){N -ki~k2 + 1){N -ki- fca)}'^' , 
dk„k2 = {fci(l + Ai)fc2(l + fc2)}'/' . 



We notice that the coefficients are defined only for < 
ki + k2 < N, which is the lower left triangular part of the 
corresponding matrix representation and the coefficient 
bkiM '^ot symmetric with respect to the exchange of 
the indexes. 

The nonlinearity, though being responsible for the very 
existence of the intraband tunneling, only defines the 
time scale and it is the lattice parameter A which en- 
ters the Hamiltonian in Eq. ( [T2| (A is a ratio of two 
integrals of the Bloch waves which are defined solely by 
the lattice). 

We conclude this section with the estimate for the en- 
ergy range: 



l-h2A 

12 
1 - 2A 



A 

2N 
A 

2iV 



<{H)< 
< (H) < 



1 + 2A 

4 

1 + 2A 



A 

2iV' 
A 

2iV' 



(16) 



(see Appendix [A| important for the numerical simula- 



tions. Since the Hamiltonian ( 12 1 is bounded from above 
and from below it follows that the energy spacing for the 
quantum particle satisfying equation (111 is, for the most 
of the spectrum, on the order of SE ^ N~'^: for a given 
number of atoms N the dimension of the Hilbert space 
is {N + 1){N + 2)/2 i.e. « in the limit N < 1 

(in the 2D case, considered in Ref. [T], the dimension of 
the respective Hilbert space was -t- 1 and respectively 
the energy distance between adjacent energy levels was 
determined by the factor N~^). 

III. THE SEMI-CLASSICAL APPROXIMATION 

A. The governing dynamical model 

The semi-classical approach employed here is simi- 
lar to that of Ref. [J. We define h = 2/N, = 
ki^2/N. Then, assuming existence of a regular function 
i^{xi,X2) 



^^^Cki,k2J Eq- (|15|) is cast as 



ihdr^j= ^a{xi,X2)^J+ ^jfoh (^xi -~'^,X2^e +hh + e'^' + h (^^2 - '^,xi ] e 

+bh (x2 + ^^^ij e'^'' + dh (xi - ^,X2 



+ ^) e'^P^'P^^ + d, (^xi + J,X2 - J) e^^P^-P^^y, 









Xi 


xi ^ 


-1) 








Xi 









(1 - Xi - X2) 1 - Xi - 2)2 + 



1/2 



where pj = —ihdx^ and we have introduced the functions, defined on the triangular domain Q < xi + X2 < 

a{xi,X2) = 1 + 2(2A - 1) [60(2^1, 3^2) + bf^{x2,xi) + d(3{xi,X2)] , 

bhixi,X2) ' ' 

dh{xi,X2) = Xi 1 Xi -t- - I 2:2 1 X2 



1/2 



(17) 
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We have = —ihSj^k, i-e. the usual canonical com- 

mutator of the momenta and coordinates. Evidently, the 
role played by h is of the effective Planck constant. The 
semi-classical dynamics corresponds to the limit ft. — > 0, 
i.e. when the number of BEC atoms N — > oo. It is im- 
portant to recall that the characteristic time t of the evo- 
lution scales as {xiN)~'^t, hence the quantity Xi-^ must 
be kept fixed. These two conditions taken together con- 
stitute the usual mean- field limit for BEC. It is important 
to mention here that the coefficients xi,2 stay bounded, 
as it follows form the definition ([s]) and the normaliza- 
tion of the wave function, implying / It^jPliySfeP ^ 1/V. 
The quantity Xi-^ ~ gN/V giving the scale of the non- 
linear time Tnl — '^t^/^nl (see Eq. (fTo|) is constant 



also in the thermodynamic limit defined as N 
constant density. 



cx) at a 



The limit /i ^ 0, if it exists, corresponds to the contin- 
uous limit of the discrete equation (17) (in this respect, 
it is similar to the WKB approach used for the discrete 
three-term relation, see Ref. lOj for details). 

In order to derive the classical equation corresponding 
to the limit /i -> we set V(a;i,a;2,r) = e»'5(^i'^2,r,/i)/'i 
for a complex action S{xi,X2,T,h) viewed as a se- 
ries S = S^°^ + hS'^^^ + 0{h^). Assuming the ac- 
tion S^^\xx,X2,t) be differentiable function we get 
the Hamilton- Jacobi equation for the classical action 
5(-')(a;i,X2,T) = 5(")(xi,a;2,T)-T/2: 



= ha{x^,X2) 
+ do{xi,X2) 



A COS 
A cos 



2A - 1 



+ 5o(a;2,xi) 



A cos 



2A - 1 



(^(f ) - ^i?) + 2A - l] ^ 7^ (^if \ ^i' 



(18) 



where we have introduced the classical Hamiltonian Ti. 
The quasi-classical dynamics is sometimes more conve- 
niently described in terms of variables Zj = 1 — 2xj and 

(j^j = Si'f'^ = p, where Xj and pj are the classical lim- 
its of the corresponding quantum variables. The Poisson 
brackets of the respective classical variables read 



numbers: b 



{4>j,Zk} = lim -[pj, 1 - 2xk] = -2(5j,fc, 

{(j)j,(l)k} = {Zj.Zk} = 0. 



(19) 



The classical variables can be associated with the quan- 
tum averages by the following correspondence: 

N N~ki 

fcl=0 fe2=0 

01 = arg((6l)2fe2) 

{JV AT-fci ~j 
E E Cl^M^kM.CkM. . (21) 
fcl=0 fc2=0 J 



difference 



i\ ^ %/7V(6(^'))* and h.j ^ The phase 

/)j in Eqs. (|2l])-([22]) is not defined if 
1 (the function^^rg" in Eq. (21 1 or 



0, 



Eq. (211 or (22 1 is 



applied to zero: Ck^^k'^ = for k'j > 1). The two pEases 
are not defined also for zi + Z2 = 0, i.e. (^'3^3) = 0, 
since in this case 63 j^*) = 0. In these cases the phases 
can be determined by taking the averages of the boson 
operators corresponding to non-zero average populations, 
i.e. in the semi-classical limit instead of {b*bk)^ one just 
takes the phase of the squared nonzero amplitude or 
{b*)\ 

The classical Hamiltonian, recovered from the 



Hamilton-Jacobi equation (18), reads 



H = -(1 - zi)(zi + z2)(Acos0i + 2A - 1) 

+ ^(1 - Z2){zi + Z2)(AcOS(/)2 -f 2A - 1) 

+ ^(1 - 2i)(l - z2)(Acos(0i - 02) + 2A - 1).(23) 



02 - ^rgiiblrbj) 

{N N-ki -\ 
E E cLkM+iMCk,M+2\-m) 
ki=0 k2=0 ) 

The first equalities in these formulae can be most easily 
established by replacement of the boson operators by c- 



As a result of Eq. ( [T9[ ) the mean-field equations of motion 
acquire the form (j — 1,2) 



dzj 2 dli. d(j)j 2 dH 
dr d(/)j ' dr dzj 



(24) 



Explicitly they read 
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zi 

Z2 



^(1 - zi)[{zi + Z2) sin(/ii + (1 - Z2)sin((/)i - (/)2)], 

^(1 - Z2) [{zi + 2:2) sin 02 + (1 - 2i)sin((/)2 - (/-i)] , 

^(1 - 2zi - z2)(Acos0i + 2A - 1) + ^(1 - z2)[cos(?!)2 - cos((?!)i - 



^(1 - 2Z2 - Zl)(AcOS02 



- 2A - 1) + ^(1 - zi)[cos(?!)i - cos((?!)i - (1)2)]. 



(25a) 
(25b) 
(25c) 
(25d) 



B. Stationary points 

Either from the point of view of dynamics governed 
by Eq. ( 25 1 or from the point of view of practical ap- 



plications the most relevant first step in studying the 
mean-field dynamics is the investigation of stationary 
points Pj = {z^x^^ .z^^^ ,<^S^^^^ and of their stabil- 
ity. We emphasize that now the stability is understood 
in the classical mechanics sense, unlike the modulational 
instability of the Bloch states mentioned in the Intro- 
duction and resulting in developing of the spatial struc- 
tures [31115]. 

We will use the notations C,a = Za — ^ and — 

(st) 

(f>a—4>a j fo'^ local coordinates describing small deviations 
from the stationary solutions. The relation between the 
populations Xj and dynamical variables z reads 



1 



(a;i,a;2,a;3) = ;r(l - zi,l- Z2,zi + Z2). 



(26) 



It is convenient to separate internal stationary points, 
i.e. the ones for which all three X-points arc popu- 
lated from the boundary stationary points for which ci- 
ther one or two X-points have zero population. The 
boundary stationary points correspond to the effectively 
low-dimensional (2D or ID) distributions of atoms in the 
3D lattice. 



1. Internal stationary points 

1. The first internal stationary point is given by 
Pi = {|,|,0,0} and describes equally populated X- 
points with zero phases (phase differences). The Hamil- 
tonian in the vicinity of Pi reads 

Hp, = i(l-3A)(C? + C| + CiC2) 



1 



(27) 



It can be diagonalized using the generation function P2 = 
(pi +P2)'fii — 2p2952, which results in (here and in similar 
formulas below we drop nonessential constant terms) 



1 



Hp, = -(l~3A)pf + -(l-3A)p^ 



A 
12 



2 
9i 



A 
36 



92', (28) 



where {pi,qi) and (p2j'?2) are the local canonical vari- 
ables. Thus Pi is unstable for A < Ai = 1/3 and corre- 
sponds to saddle points on the planes {pj,qj), while it is 
linearly stable otherwise, though corresponds to a local 
maximum of the Hamiltonian. The respective motion 
can be interpreted as a 2D linear oscillator with nega- 
tive effective masses and with equal frequencies f2p, = 
[A(3A- 1)/12]5. 

We observe that the critical value of the lattice parame- 
ter Ai for equally populated X-points coincides with that 
in the 2D optical lattice (see [UIJ). 

2. The second stationary point Pj = {|'|'^'~^} 
also corresponds to equal populations of the AT-points, 
but characterized by mutual 27r/3-phase differences. In 
this case the Hamiltonian about the stationary point 
reads 



1 



A 



Hp, = ^{2-SA){Ct + C2+CiC2) + ^ivl + ^l 



18 



A 



[Ci(2V2 



A 



Vl) + C2('y52 - 2^1)] + ^ - 3 



^1^2) 
1 



The variables and ip are now mixed and the transforma- 
tion which diagonalizes the Hamiltonian is complicated. 
The eigenfrequencies, however, can be directly obtained 
by considering the characteristic equation. We get two 



distinct values il 



(±) 



± 



6 



4 



which become 



complex for A > A2 = 2/3. Therefore, the P2-point is 
linearly stable for A < A2 and is unstable otherwise. The 
critical value A2 is a characteristic feature of the 3D case 
and does not exist in the 2D setup. One also readily con- 
cludes from the symmetry that another stationary points 
isgivenby P^ = {i,i,-f , f }. 

3. The next three stationary points are given by 

= {'A'A-^n}, P4 = {l±A,3A5^,0,7r} and 

P^ = i±A^7r,o|. The point P3 corresponds to 

the populations 



{xi,X2,X3) 



1 -A 1 - A 1 + A 



3-A'3-A'3-A 



while the points P4 and P4 correspond to the 
same distributions with X-points being interchanged: 

{xi,X2,X3)p^ ^ f^'i^i'l^) and (a;i,a;2,X3)p' = 
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(l^S' l^)' ^^'^ stability properties and the diag- 
onalized local Hamiltonian are the same for these three 
points. Consider, for instance, P3. The local Hamilto- 



equally populated (see Appendix B). It reads: 
'1 1 



{xi,X2,X3) = 



2' 2 







(33) 



+ 4(1 - m! + CD + 4(1 + A)CiC2 



A(l-A) 
(3-A)2 



[A(^? + </5^) + (l-A)^iH (29) 



can be diagonalized by means of the canonical transfor- 
mation generated by F2 = {pi + P2)vi + {Pi ~ P2)^2'- 



np. = J(3 - A)p? + 1(1 - 3A)p^ + "^-^ ql 



A(1-A)(1-3A) 2 



4(3 - A)2 



92- (30) 



Hence P3 is a saddle point of the Hamiltonian in the 
plane (j32,92), and thus is always unstable (the same is 
true for P4 and P4). 

4. In the critical case, A = Ai, for the zero phases 



^'1,2 = the classical Hamiltonian (23 1 is flat in (zi, Z2) 



7i(zi, Z2, 0, 0) = and zi_2 — due to the phases, i.e. 
the whole domain of (zi, Z2) has the same energy for zero 
phases. 



2. Boundary stationary points 



have 



where cos{4)\) ~ ^^x"- terms of z we 

= {OiO:0A:0a}- Noticing that this point is also a 
stationary point of the Hamiltonian ( 23 1 we obtain the 
local Hamiltonian as follows 



Hi 



i(l-3A)(2 + Ci' + C2D-^('Pi-^2) 



A 



sin (/)a(Ci + C2)(<y5l + (P2) 



(34) 



Next, using the local canonical transformation with the 
generating function F2 = (Pi + P2)vi + (Pi - P2)v'2 we 
arrive at (dropping the constant) 



4(l-3A)(p^ 



-sm0APigi - 



-(1-3A) 



pI 



pI 



A^ sin 2 
(1-3A)2'^\ 



A 



<Z2^ 



(35) 



with pi = pi - 



A sin <pA „ 

1-3A yi- 



Therefore, Pb^ is a saddle point 



in the whole domain of its existence (except in the critical 
case A — 1/3), hence it is unstable unless A = Ai. 

For the sake of convenience, in Table |l] we present the 
list of the stationary point and their linear stability prop- 
erties. 



As it was mentioned above, there exist boundary sta- 
tionary points corresponding to all atoms populating only 
one X-point or only two X-points. As the phases be- 
come undefined in such a case, it is convenient to use 
the semi-classical Hamiltonian obtained directly from 



Hamiltonian 



^ by 



the substitution b 



7V6*l-'='^ and 



' s/NV^^'' (to have normalized amplitudes). 

1. Consider the solutions with b\ 2 — f3i.2 and 63 — 
I + P3, 1(3 j \ ^ 1, i.e. close to the stationary point Pg^ de- 
scribing all atoms occupying just one X-point (X3-point 
in this case): 



Using that ll-H/Jgp = 1 



(a;i, 0:2, 2:3) = (0,0,1). (31) 
|/32p we obtain the local 



coordinates 

{zi,Z2,<j}l,<j}' 




stability 



stable for A > | 
stable for A < I 
stable for A < | 
unstable 

unstable 

unstable 

stable for A < | 
unstable for A 7^ 



populations 

{X1,X2,X3) 



(1 1 1) 

\3' 3 ' 3/ 
(-1 1 1\ 
\3' 3 ' 3/ 

/I I I) 
v3 ' 3 ' 3 / 

1-A 1- 
3-A' 3- 



1+A 
3-A , 



1-A 1+A 1-A ' 
3-A ' 3-A ' 3-A , 

1+A 1-A 1-A ' 
3-A ' 3-A ' 3-A , 



(0,0,1) 

(jjo) 



TABLE I: The stationary points of the classical Hamiltonian 
and their linear stability properties 



Hamiltonian as follows 



n 
n 



X3 

(i) 



n 



(1) 
X3 



n 



(2) 
X3' 



*^2 



+ /?■]■ (32) 



IV. QUANTUM EVOLUTION 



A. The initial state and the numerical approach 



The eigenfrequencies are equal for the two modes (3iX- 
rig^ = i[3A2~4A-f 1]5. Hence for A < Ai this stationary 
point is a local minimum and is stable, while for A > Ai 
it is unstable. 

2. Moreover, it is easy to show that for A > Ai there is 
one more stationary point with two X-points being 



We are interested in quantum dynamics of an ini- 
tial state of a large number of atoms with the aver- 
age values of the occupation numbers and the phases 
being close to a semi-classical stationary point. To 
find out the structure of such an initial state con- 
sider the semi-classical wave function 'ip{xi,X2,T) = 



8 



^iS(xi,x2,T,h)/h ^ where in the lowest-order approximation: 
S{xi,X2,T,h) = S'^'^^^xi, X2, t) + 0{h) . In the vicinity of 

a stationary point {x'f''\x'^^^) we can expand the classi- 
cal action as follows 



„(c^)^2 



1 

+O[(xi-4^'0' + (2^2-4"01- (36) 



"2 

„(c^)^2l 



Therefore, recalling that Xj — kj/N and h — 2/N and 
taking into account that the average values of Xi^2 niust 

be close to the semi-classical ones x^i2i can approxi- 
mate the initial state by the Gaussian function 



ki ,k2 



Coe' 



i)\ 'fel+02 '=2)-^ '-—T — 



.(37) 



Here ^^^^2 ^^'^ ^'^'^ ^^'^ classical phases and popula- 
tions, Co is the normalization factor and (Tat is the width 
parameter such that 



(38) 



(the first inequality is imposed to guarantee smoothness 
of S{xi,X2, T, h) with respect to 2:1^2 and the second one 
is the condition of small width of the wave-packet in the 
Fock space). Due to symmetry of the quantum Hamil- 
tonian (bosons are created by pairs), the classical phases 
01,2 give rise to six different quantum states of the form 
(37l with the phases <pi^2 + 27rsi_2, si,2 G {—1, 0, 1} (see 



also the discussion of phase states below) . 

One can expect that the state ( 37 1 , ( 38 1 with the classi- 



cal variables satisfying the respective Hamiltonian equa- 
tions is a good approximation for the actual quantum 
state for all times r as /i ^ if the classical stationary 
point is stable. Indeed, in this case the expansion (36) 
can be truncated as indicated. 

To get a numerical solution of Schrodinger equation 
with a controllable accuracy we have used the 
method of Ref. i.e. the expansion of the unitary 

operator U = exp{—iNHT} over the Chebyshev polyno- 
mials 



(151 



-iNHAr 



-iNEt 



Ce{NAEAT)Te{I), (39) 



where E = (E^,^^ + E^i,,)/2, AE = (S^^ax - 
£^„iin)/2, with E'min and i?max bcing the lower and up- 
per bounds taken from equation (16 1, Ti{I) being the 
£-order Chebyshev polynomial of the Hermitian opera- 
tor I = {H — E)/AE with the eigenvalues lying on the 
interval [—1,1]. The coefficients are given as Ci{>c) = 
(— i)^(2 — Se,o)Je{>c) where Je{>c) is the Bessel function 
of the first kind. Due to the uniform convergence of the 
Chebyshev series on [—1, 1] and the fact that the coeffi- 
cients vanish exponentially for sufficiently large £ (for a 
fixed At) one can compute the evolution operator for the 
Schrodinger equation at the times r — At, 2At, 3At, . . . 
with arbitrary given accuracy, limited only by the round- 
off errors (we have set the error to be of the order 10~^). 



B. Quantum evolution about the Pi-state 

For A < Ai Pi is a saddle point and is unstable with re- 
spect to small perturbations. An initial quantum state in 
the form (|37|) such that the average initial populations Xj 
and the phases cj)j are close to the semi-classical station- 
ary values Xj = 1/3 and (j)j = results in the evolution 
presented in Fig. [T] The initial localized, nearly-Fock, 
state transforms to a broad oscillating state (lower pan- 
els of Fig. [1]) persisting at least for some long evolution 
time. 
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FIG. 1: Quantum evolution of A'' = 200 BEG atoms loaded 
into the high symmetry points Xi, X2 and X3. The lat- 
tice co nst ant is A = 0.21. We use the initial state as 
in Eq.(37l with a = ^ N/2 with the initial populations 
{xi,X2) = (0.330,0.337) and phases {<j)i,(t>2) = (0.02,-0.02). 
The initial stage of evolution is given in the upper two panels, 
while the lower two panels show an oscillating state by which 
the initial (localized) state is replaced. 



The emergent state can be approximated by a linear 
combination of a small number of the phase states. The 
one-dimensional phase states are defined here via the dis- 
crete Fourier transform (DFT) 



1 ^ 



(40) 



where 9e = Evidently {6(/\9e) = Se>^i. Therefore, 

the phase states give another basis of the Hilbert space, 
in fact 



(41) 



For a fixed total number of atoms rii + 712 + 113 = N 
the three-dimensional phase states are projected onto a 



9 



two-dimensional subspace, i.e. the wave function can be 
written as 



N N 



(42) 



where 



=0 i2=0 



N N-ki 



N 



^ E E e-^'^^'i'^'^'^^Cfc^fc, (43) 



fei=0 k2=0 



is nothing but the DFT of the coefficients Ck^^k^ extended 
over whole domain of < fci^2 < by padding them with 
zeros. Note that the phase 9 is half of the value of the 
semi-classical phase </> in the limit ft, — > 0. 
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FIG. 2: The DFT transform of the wave function of Fig. [T] 
at two large times (two panels are used to show the relatively 
small deformation with time). 



We find that the DFT of the wave function of Fig. 
[l] is concentrated at the following values of the phases 
9i = {0,±7r}, see Fig. |2] These states correspond to 
the phases in the semi-classical limit, i.e. to the 

phases of the stationary point Pi. 

The quantum evolution of the initial state correspond- 
ing to a stable classical stationary point as ft is dif- 
ferent, see Fig. [S] First of all, the localized (i.e. nearly 
Fock) state remains localized. Note that the quantum 
and the semi-classical dynamics are very close in this 
case, see Fig. |4] though the number of atoms is rather 
small. 

In the nonlinear tunneling of BEC in a square 2D op- 
tical lattice yy (where the two-mode model appears) the 
quantum evolution features appear as collapses and re- 
vivals of the semi-classical dynamics. The energy spacing 
5E ~ iV~^ discussed in in Sec. |ll] for the three-mode 
model (as compared to 5E ^ N^^ for the two- mode 
model) prevents observation of the quantum collapse. In- 
deed, the semi-classical regime requires large number of 
atoms, thus large evolution times r ~ iV^ are required 
for observation of the first quantum collapse. We verified 
that the quantum oscillations of Fig. |4] follow the semi- 
classical ones without occurrence of the quantum collapse 
for times up to r = 60000 at least, which would exceed 
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FIG. 3: Quantum evolution of an initial Gaussian state of 
TV = 200 BEC atoms with a = \/iV, the initial populations 
{xi,X2) = (0.35,0.32), and phases ((^1,(^2) = (0.01,-0.02). 
The lattice constant A = 0.36, i.e. the semi-classical state Pi 
is stable. Oscillations of the wave function about the initial 
state are observed (the time increases clock- wise). 
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FIG. 4: Comparison of the quantum evolution of Fig. |3]with 
the semi-classical evolution corresponding to the initial av- 
erage values of the populations and phases. The upper and 
lower panels show average populations and phases, respec- 
tively. 



by far the lifetime of BEC (see Sec. [IT]). This result also 
suggests that the quantum collapse may not exist in the 
model at all. 



C. Quantum evolution about the P2-state 

The above results show that the quantum model of N 
identical bosons distinguishes between the stable and un- 



10 



stable classical stationary points. This conclusion agrees 
with the correspondence between the quantum stability 
of a semi-classical state in a system of identical bosons 
and the Hamiltonian stability of the corresponding sta- 
tionary point in the classical limit il21 . 
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FIG. 6: Comparison of the quantum evolution of Fig. [s] (solid 
lines) with the semi-classical result corresponding to the ini- 
tial average values of the populations and phases. The upper 
panel gives the average populations and the lower one the 
phases (dashed lines). Top panel gives the population x\ and 
the bottom one the phase <^i. 
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FIG. 5: Recurrence of the wave function in quantum evo- 
lution (the time increases clock-wise) of an initial Gaussian 
state of A'' = 200 BEG atoms with cr = ViV, the initial 
populations {xi,X2) = (0.330,0.337), and phases (0i,(^2) = 
7r(2/3 + 0.01, -2/3 -|- 0.006). The lattice constant A = 0.64, 
i.e. the semi-classical state P2 is stable. 



However, due the discreteness of the quantum energy 
levels the quantum evolution can have features not found 
in the classical model (the two cases, of course, agree in 
the limit N 00 when the quantum energy spacing 
goes to zero). This is clearly illustrated by the results 
presented in Figs. |5] and |6] Indeed, Fig. [5] illustrates 
one period of the wave-function spread and subsequent 
recurrence to the localized distribution, which is respon- 
sible for the deviation of the quantum averages from the 
corresponding classical variables, see Fig. |6] Note how- 
ever, that the quantum averages remain close to the the 
classical stationary point values, in accordance with the 
general correspondence of the quantum and classical sta- 
bility [15]. 

One more difference is apparent in Fig. [6] as compared 
with Fig. |4] the semi-classical dynamics about the P2- 
point features two frequencies instead of one, as it is for 
the Pi-point. Despite the disagreement of the quantum 
averages and the classical dynamics, the recurrence pe- 
riod is in fact very close to one of the classical oscillations 
periods t 50. 

The quantum dynamics corresponding to the unstable 
classical fixed point P2 is similar to that in the case of 
unstable Pi-point, namely the localized, i.e. nearly Fock- 
state, is replaced by a linear combination of a small frac- 
tion of the phase states, see Figs. [7] and [Sj The phase 
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FIG. 7: Quantum evolution of an initial Gaussian state of 
TV = 200 BEG atoms with a = vTS, the initial popula- 
tions (xi,X2) — (0.3300.337), and phases {4>\,4>2) ~ 7r(2/3 -I- 
0.07,-2/3 -I- 0.07). The lattice constant A = 0.69 and the 
semi-classical state P2 is unstable. 



states of Fig. [8] are concentrated about the following 
phases: 

J / 27r tt\ ( 2-K 27r\ / tt ttn /tt 27r\ 

which correspond to the classical phases (0i,(/)2) = 
{(±^,=F^)}, i.e. to the phases of the stationary point 
P2 and its equivalent Pj- 
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FIG. 8: The DFT transform of the wave function represented 
in the right panel of Fig. [7] 
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For a special initial atomic distributions it is possible 
to have stable-like quantum dynamics about an unsta- 
ble semi-classical fixed point which is conditionally stable 
for the special initial conditions. For instance, the fixed 
point P3 is stable for the initial states with no (^2,92)- 
components in the classical limit (which is supposed to 
be a small perturbation about the fixed point). In this 
case the wave function remains localized and performs 
oscillations about the initial state (not shown). 

As the stationary points corresponding to unequal pop- 
ulations of the X-points of the lattice are unstable and 
loading BEC into the unequal distribution among the 
high-symmetry points is not an easy (if at all possi- 
ble) task we discard the further analysis of the dynamics 
about the points P3, P4 and P4. 



D. Dynamics of the boundary states 

One can easily load BEC into a single X-point by 
switching on a moving cubic lattice. Thus, it is important 
to consider the boundary stationary point Pbi ■ Let us 
consider Xs-point being initially populated. For A < 1/3 
the point Pb^ is classically stable and the quantum dy- 
namics consists of localized oscillations about the initial 
state. If however, the lattice parameter passes the critical 
value Ai the instability of Pg^ results in tunneling to the 
equal distribution of atoms between the three X-points, 
see Fig. [9j 

The dynamical instability of the Pbi can be used to 
prepare the system in the equal distribution of atoms 
between the A"-points by loading first the Pg^ -point as 
discussed above and modifying the lattice parameter A to 
force the dynamical instability of Pg^ to develop, i.e. as 
shown in Fig. [9] The oscillations in Fig. |9]are about an 
equal distribution of atoms between the X-points and 
the zero phases, thus the quantum state is the semi- 
classical state about the Pi-point of the form given by 
equation ( [37| (or a linear combination of such states). 
Moreover, one can notice that the energy of a stationary 
semiclassical state {N ^ 1) in the main order is given 



FIG. 9: Comparison of the quantum evolution (solid lines) of 
with the semi-classical one (dashed lines). The upper panel 
gives the average population xi and the lower one the phase 
01. Here the initial populations and phases are {xi,X2) = 
(0.64, 0.04) and (<^i, ^2) = (0.05, 6), the lattice parameter A = 
0.41, N = 260, and a = \/N. The stationary point Pb^ 
(xa — 1) is unstable. 



by the zero-point energy of the local classical Hamilto- 
nian, since the energy spacing between the local bound 
states is on the order or smaller than 1/7V (since the 
quantum oscillator model, obtained by the "reverse quan- 
tization" procedure of the local classical Hamiltonian, 
has the energy spacing 0(h)). Comparing the energies 
Ep^^p, = A/2-l/3-fO(l/7V) and Ep, = -l/3+0{l/N), 
we see that the ground state for A > 1/3 corresponds to 
the Pi-point. 

On the other hand, the P2-point and its equivalent 
point P2 are stable and Pi is unstable for A < 1/3 
(see table |l|. Note also that the zero-point energy of 
Ep^ = A/2 < 1/6 is lower than that of Ep^^ = 1/4, an- 
other stable point for A < 1/3. Thus, given the quantum 
state with an equal distribution between the AT-points, 
i.e. Pi, one can prepare another such stable state (in fact 
P2 or P2 or their linear combination) by repeating the 
above procedure but now starting from the Pi-point by 
adiabatically changing the lattice parameter to A < 1/3 
followed by the thermal cooling procedure. 

Stationary point Pg^ is unstable in the domain of its 
existence A > 1/3 (except for the critical value A = 1/3). 
This stationary point corresponds to the quasi 2D sta- 
tionary state, however its instability rules out observation 
of 2D quantum dynamics P], for instance the quantum 
collapses and revivals. We have found that an initial state 
with almost equal distribution of atoms between two Ap- 
points results in the sequence of quantum recurrences, 
when the wave function returns to a state with almost all 
atoms distributed among the initially populated points, 
see Fig. 10 (in the figure this state corresponds to an 

X2 = l). 
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FIG. 10: Quantum evolution of a state corresponding to ini- 
tially almost equally populated two X-points (here Xi and 
X2). The time increases clock- wise. Here the initial popu- 
lations and phases are {xi,X2) = (0.48,0.48) and (01, 02) = 
(;!>a(0.98, 1.03), the lattice parameter A = 0.51, N = 200, and 
a = \/iV. 



V. DISCUSSION 



lowed by the quantum recurrence to the initial nearly 
Fock state. This is reflected in a deviation of the quan- 
tum averages from the semi-classical dynamics. In the 
case of an unstable stationary point, the initially local- 
ized state, i.e. nearly Fock state, is replaced by a nearly 
phase state with the phases concentrated at the semi- 
classical value corresponding to the unstable point (more 
precisely, a linear combination of nearly phase states, 
since the quantum phase appearing in the wave function 
of the effective quantum particle is equal to half of the 
semi-classical phase due to the symmetry of the quantum 
Hamiltonian) . 

Existence of the stable stationary point with all atoms 
populating just one X-point of the lattice allows for the 
experimental study of the 3D nonlinear tunneling by 
modifying the optical lattice to change the value of the 
lattice parameter A. When the instability of the singly- 
populated X-point is reached by varying A, the quantum 
evolution quickly establishes equal distribution between 
the three degenerate X-points. 
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The nonlinear tunneling of BEC with a large num- 
ber of atoms N can be considered in the semi-classical 
approximation, with the effective Planck constant being 
h = 1/N. We have considered the correspondence be- 
tween the semi-classical regime (equivalent to the mean- 
field regime) and the full quantum regime of nonlinear 
tunneling between the AT-points of the Brillouin zone 
of a cubic 3D lattice. In particular, we have derived a 
quantum three-mode model and rewritten it as a two- 
dimensional Schrodinger equation for an effective quan- 
tum particle, where the effective Plank constant is 1/A^, 
the time scale is determined solely by the nonlinearity 
of BEC, while the dynamics is controlled by a lattice 
parameter A. The corresponding semi-classical model 
is the mean-field approach taking into account the oc- 
cupations of the AT-points only. Though we have used 
rather small number of atoms, N = 200, we have found 
the regimes of excellent correspondence, these are mainly 
about the stable stationary points of the mean-field ap- 
proach. In particular, numerical simulations show that 
the quantum dynamics about the semi-classical station- 
ary point distinguishes the stable and unstable cases. In 
the case of a stable semi-classical point, one scenario con- 
sists of the wave function performing oscillations about 
the initial state with the averages following the semi- 
classical dynamics. The discreteness of the quantum en- 
ergy space, however, leads to a scenario not present in the 
semi-classical case: the sequence consisting of the wave 
function spread (i.e. becoming a nearly phase-state) fol- 



APPENDIX A: THE ESTIMATE ( 16 ) 



We will use the inequality 

m]rbl + {blfb'^)\ < {n,{n, - 1) + n,{n, - 1)), (Al) 
which follows from 

< {^\{A + B)\A + B)\^) 



by setting A = bj and B = ±6^. Using ( Al I we obtain 



{H-) < {H) < {H+), (A2) 
where H± are two c-number operators: 



rij Uk A X ^ Uj Uj — 1 
iVlv 2"^ iV N 



j = l ' ■ j<k " - - - j = i 



which may be treated as classical functions of Uj. Next, 

reducing the total squares (X]^=i ^j/-^)'^ = 1 in H± we 
have 



i7_ 



4 ^ 7V2 2 
1 - 4A A A 



4 2 



N 
1 

N 



(A3) 
(A4) 
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The inequalities (16 1 follow from (|A3| and (A4) if one By changing cos 9 sin 6* and (f>i 



we have also 



takes into account that N^/3 < J2j=i "-j 



APPENDIX B: THE SEMI-CLASSICAL 
STATIONARY SOLUTIONS WITH ONLY TWO 
POPULATED X-POINTS 



Consider the Hamiltonian ( 12 1 in the semi-classical 
limit, i.e. bj — > \fNh''^^\ Suppose that there is a sta- 
tionary point with 6io — cos ^e""^^/^, 620 — sin 6'e~*'^^/^, 
and 630 = with some Q < 9 < tt /2. To find all such 
stationary points we use that the semi-classical Hamil- 
tonian expanded about such a point does not have any 
linear terms. Without loss of generality we can set 

hi = cos9e-"l'^^^{\ + l3i),h2^siTi9e-'"l'''^{\ + 132), 
&3 - /33 e Re (Bl) 

where, due to the conservation of the number of atoms, 
in the linear order we get 

Pl^- cos^ 9{Pi + (31) - sin^ 9{p2 + ^2)- (B2) 

Using this we obtain the linear in [3i 2 terms as follows: 

"cos2 6l 



dn 

Wi 

A 



cos^ 9 



Pl,2=0 



A cos^ 9 



-I sm^^e^f"^^-"^!' 



cos^ 9 cos 01 — sin^ 9 cos 02 



(B3) 



In calculation of ( B3 1 we have used that the contributing 
terms are 

Hfi^ = ^nf + A[nin2 + nin^ + n2n3] 

+ jm)X + ib*s)X + {bi)X + ib;)'bi]. 



dn 

W2 

A 



= sm a 

/3i.2=0 

cos^ 6'e'("^'2-'^'i) 



sm 



Asin^6' 



sin^ 9 cos 02 ~ cos^ 9 cos 0i 



(B4) 



The r.h.s.'s in equations (B3), (B4| and in their com 



plex conjugates should give zero for a stationary solution. 
First of all, from (B3) we have sin(02 — 0i) = 0, hence 
02 = 01 or 02 = 01 ± TT. Then, combining equations 



(B3), (B4l we get 



(cos^ 9 



1 A ^ / 
:^ ^ A — — COS! 

2 2 ^ 



) 



i.e. (i) cos 9 = sin 9=1/ \/2 or (ii) cos 
In case (i) we obtain the phase cos((/ 
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while 

in case (ii) there is no solution except for the special 
value of the lattice parameter A = 1/3 when the phase 
01 becomes arbitrary. Hence, we arrive at the stationary 
point {xi,X2,X3) = (|, |, 0), 01 = 02 = 0A7 which exists 
only for A > Ai = 1/3 (the phases appear in the initial 
state (37l). In terms of z- variables the stationary point 
reads zi = Z2 = 0. 
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